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ABSTRACT 

In N-body simulations the force calculated between particles representing a given 
mass distribution is usually softened, to diminish the effect of graininess. In this paper 
we study the effect of such a smoothing, with the aim of finding an optimal value of 
the softening parameter. As already shown by Merritt (1996), for too small a softening 
the estimates of the forces will be too noisy, while for too large a softening the force 
estimates are systematically misrepresented. In between there is an optimal softening, 
for which the forces in the configuration approach best the true forces. The value of 
this optimal softening depends both on the mass distribution and on the number of 
particles used to represent it. For higher number of particles the optimal softening 
is smaller. More concentrated mass distributions necessitate smaller softening, but 
the softened forces are never as good an approximation of the true forces as for not 
centrally concentrated configurations. We give good estimates of the optimal softening 
for homogeneous spheres, Plummer spheres, and Dehnen spheres. We also give a rough 
estimate of this quantity for other mass distributions, based on the harmonic mean 
distance to the kth neighbour (k — 1, .., 12), the mean being taken over all particles 
in the configuration. Comparing homogeneous Ferrers ellipsoids of different shapes we 
show that the axial ratios do not influence the value of the optimal softening. Finally 
we compare two different types of softening, a spline softening (Hcrnquist & Katz 
1989) and a generalisation of the standard Plummer softening to higher values of the 
exponent. We find that the spline softening fares roughly as well as the higher powers 
of the power-law softening and both give a better representation of the forces than 
the standard Plummer softening. 

Key words: galaxies: structure - galaxies: kinematics and dynamics - methods: 
numerical. 



1 INTRODUCTION 

iV-body codes are often used to simulate the dynam- 
ical evolution of galaxies and galaxy systems, even though 
the number of particles that present day computer hardware 
and software can handle is smaller, by several orders of mag- 
nitude, than the number of stars in even a small galaxy. 
Because of this the particles do not represent individual 
stars, but should be considered as Monte-Carlo realisations 
of the mass distribution in a galaxy. In such simulations, 
when close encounters between individual particles are of 
no relevance to the physical problem under consideration, 
the gravitational force between two particles is smoothed, 
in order to reduce the spurious two-body relaxation due to 
a number of particles necessarily much smaller than the to- 
tal number of stars in the system. Although a large softening 



will ensure a low relaxation rate, we can not increase this 
value at will, since a high value of the softening introduces 
other drawbacks and in particular a bias to the gravitational 
force. 

The question we will address in this paper is what value 
of the softening should be used in order for N particles to 
represent best a given density distribution. In a recent paper 
Merritt (1996, hereafter M96) argues that too small a value 
for the smoothing will give too noisy estimates, while too 
large a value will give a systematic misrepresentation of the 
force. M96 and Athanassoula el al. (1998; hereafter A+98) 
found that this optimal value of the softening depends on 
the number of particles N approximately as N~ 0,3 . In this 
paper we will extend previous work to other configurations, 
other Af-body algorithms for calculating the force and other 
functional forms for the force approximation. The two N- 
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body methods under consideration are briefly described in 
section [2], where we present also the tools we will use for our 
comparisons. The optimal softening for the case of a Plum- 
mer sphere is discussed in section H, including comparisons 
of MISE and MASE, as well as comparisons of calcula- 
tions with GRAPE-3 and GRAPE-4. Section [5| presents the 
case of two Plummer spheres and section ^ two other density 
distributions, one more and the other less concentrated than 
the Plummer sphere. Section ^ discusses the effect of triaxi- 
ality, using the Ferrers' ellipsoid. Different types of softening 
are introduced in section [?] and a different way of calculating 
the force, using a hierarchical octal tree, in section ^. Finally 
our results are summarised and discussed in section H. 



2 METHODS 

2.1 Force evaluation methods 

Several methods have been developed for calculating the 
forces in self-consistent iV-body simulations (see e.g. the 
reviews by Sellwood 1987 and Athanassoula 1993). They 
all have their advantages and disadvantages and which one 
should be chosen depends to a large extent on the problem 
to be solved. We will here briefly describe two methods not 
using either grids or expansions, in which case the introduc- 
tion of the softening is more straightforward. 

1) Direct summation. This method involves calculating 
the forces between every pair of particles and then summing 
up all contributions to the force on a given particle. It is 
straightforward, easy to program and easy to vectorise and 
parallelise. For a fixed number of particles it is by far the 
most accurate method, since its only approximation is the 
introduction of softening, and can be used without any re- 
strictions on geometry. This method was used heavily in the 
seventies and early eighties, but was soon considered as a 
dead end because of its large claims in CPU time. Indeed 
the amount of time necessary for one force calculation scales 
as N(N — 1), where N the number of particles in the simu- 
lation. Thus, until recently, direct summation could not be 
used freely for simulations where a sufficiently high number 
of particles is necessary, as e.g. simulations of disc galaxies. 
Two major advances in computer hardware have changed 
the situation in the past few years. 

The first is the advent of parallel machines, either SIMD 
or MIMD, which, given the simplicity of the communications 
involved in the direct summation method, make it possible 
to achieve relatively high performances with little software 
investment. A large number of such systems are actually 
in use, from the powerful many-node CRAY T3E, to "Be- 
owulf" clusters, which have been recently implemented in 
many universities and research centers. 

The second is the realisation of GRAPE, a dedicated 
computer card which performs the force calculation by di- 
rect summation and which can be coupled to a standard 
workstation allowing one to achieve at relatively low cost 
an excellent CPU performance. A series of such GRAPE 
boards have been built by the group in Tokyo University, 
starting with GRAPE-1 and evolving steadily to GRAPE-5, 
while new members of this family are under development. 
For a brief history of this project see Makino & Taiji (1998) 
and references therein. Boards with even numbers have high 



accuracy arithmetic and can be used for collisional simula- 
tions, where close encounters play an important role in the 
evolution of the system, as for globular clusters and planetes- 
imals. Boards with odd numbers have a more limited pre- 
cision arithmetic and can only be used for collisionless sys- 
tems. Thus GRAPE-3 uses 14 bits to represent the masses, 
20 bits for the positions and 56 bits for the forces. Neverthe- 
less this was shown to be sufficient for simulations (A+98), 
since the dominant source of error is the noise in the form 
of two-body relaxation, while the effect of the error in the 
force calculations is comparatively smaller (Hernquist et al. 
1993, Makino 1994). 

Two such systems have been used for the calcula- 
tions presented in this paper, namely the GRAPE-3AF 
and GRAPE-4 systems in Marseille Observatory. The for- 
mer consists of five GRAPE-3AF boards coupled via an 
Sbus/VMEbus converter to a Sun Ultra 2/200, and is 
described in detail in A+98. The latter consists of one 
GRAPE-4 processor board and one control board (Makino 
et al. 1997) linked to an Alpha 500/500 workstation via a 
PCI interface board (Kawai et al. 1997). 

2) Treecodes. With the help of treecodes considerably 
more particles can be used in JV-body simulations, while one 
of the main advantages of direct summation codes, namely 
that they can be applied to systems with any geometry, 
is preserved. They stem from the simple idea that when a 
group of particles is sufficiently distant from another particle 
and its spatial extent is small with respect to the distance 
separating it from the particle, then this group can be con- 
sidered as one entity, and only monopole and, in some cases, 
quadrupole terms need be retained for the calculation of the 
force exerted by this group on the particle. Whether a group 
of particles can be considered as one entity, or whether it has 
to be further subdivided in subgroups is determined by the 
tolerance, or opening angle, which determines the precision 
of the force calculation. In this way one obtains a consider- 
able saving over the number of force calculations necessary 
in a direct summation code, so that for tree-codes the neces- 
sary time increases with the number of particles N as NlogN 
or as N. The version of the treecode most commonly used in 
astronomy is the Barnes- Hut algorithm (1986) and in par- 
ticular its vectorised version, freely available from Hernquist 
(1987). 

The treecode shares several of the advantages of di- 
rect summation, namely that it can be applied to distribu- 
tions with any geometry and large spatial and/or temporal 
density variations, while being considerably faster. It has, 
however, also many drawbacks. It is more difficult to pro- 
gram, and its vectorisation, and particularly parallelisation, 
present several difficulties. A more worrisome drawback is 
the fact that Newton's third law of motion is not necessar- 
ily obeyed. Take the example of an isolated particle A and a 
particle B in a cluster of particles far from A. Then the force 
of A on B is correctly calculated, while A will only feel the 
effect of the cluster of particles as a whole. Even so the total 
force on A from the whole of the cluster is adequately cal- 
culated and therefore one can expect a correct evolution of 
the system. Other drawbacks are discussed by Salmon and 
Warren (1994), who discussed the relative merits of opening 
angle criteria (or "Multipole Acceptability Criteria" , as they 
call them). 

The implementation of such a code on GRAPE systems 
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is possible, provided one keeps the monopole term only, and 
neglects quadrupole and higher order terms (Makino 1991, 
A+98). Such terms are nevertheless possible to include using 
a method described by Kawai & Makino (1999). 

2.2 Notations and computing miscelanea 

The usual way of softening the Newtonian force exerted 
on a particle by another particle is to by replace in the cal- 
culation of the potential the distance between the two par- 
ticles, r, by \/r 2 + e 2 , where e is the softening parameter. In 
this way the acceleration (or force per unit mass) on particle 
i from TV — 1 other particles is written as: 



Fi = G 



N 

^/(|X, -S.|- + 
3=1 



(1) 



Here G is the gravitational constant which we hereafter take 
to be equal to unity. In the following we will assume for 
simplicity that all particles have the same mass, equal to 
m. Each term in eq. (1) is equal to the force felt by a point 
of unit mass in a Plummer sphere of scale-length e, and we 
shall thus hereafter refer to this softening for brevity as the 
Plummer softening. It was initially introduced by Aarseth 
(1963) in simulations of clusters of galaxies. Although this 
is the most commonly used form of the softened Newtonian 
force, it is not the only one. Alternative types of softening 
will be described in section 

Independent of the way softening is introduced, it is 
clear that too little smoothing leads to noise and too much 
of it to the modeling of a gravity which is far from Newto- 
nian. In order to compare the calculated forces to the true 
ones M96 introduced the quantities ISE (Integrated Square 
Error), ASE (Average Square Error), MISE (Mean Inte- 
grated Square Error) and MASE (Mean Average Square 
Error). We will use them here also, since they provide a very 
useful way of quantifying how well the force of a given con- 
figuration is approximated. We will, however, modify them 
somewhat for our present needs, e.g. by introducing a multi- 
plicative constant to make them dimensionless, so that it is 
possible to make comparisons between models with different 
total mass and size. We briefly summarise all the definitions 
below. 

Let F trtle (xi) be the true force from a given mass dis- 
tribution at a point x^, and let Fi be the force calculated 
at the same position from a given TV-body realisation of the 
mass distribution and using a given softening and method. 
Then the average deviation between the two forces is given 
by 



JV 

ASE = — ^ |Fi - F true ( Xl )| 2 



(2) 



where the summation is over the TV particles in the real- 
isation. We can similarly introduce the "integrated square 
error", or 



ISE =M 



p(x)|F(x)-F true (x)| 2 dx 



(3) 



where p(x) is the true density at point x, M is the total 
mass in the system, F(x) is the force at position x calcu- 
lated from the TV-body realisation of the mass distribution, 



and the integral is taken over a volume in space encom- 
passing the configuration. For spherically symmetric mass 
distributions using the ISE rather than the ASE brings a 
substantial gain in CPU time. Indeed the integration over 
the angles can be done analytically and one is left with a 
one-dimensional integration along a radius. We will some- 
times refer to this as the radial ISE algorithm. For this the 
gain in time is proportional to ni nt /(TV — 1), where riint is the 
number of points used when calculating the integral. Nev- 
ertheless this is made at the expense of considerable noise, 
since the number of points at which the integrand is cal- 
culated isrelatively small. This will be discussed further in 



section 3.3 



In order to get rid of the dependence on the particular 
configuration, which is of no physical significance, we gener- 
ate many realisations of the same smooth model and average 
our results over them. Thus the mean value of the ASE is 



c N 

MASE = — < V 
TV ^ 



Fi — Ft rue (xi)| > 



(4) 



where <> indicates an average over many realisations. Sim- 
ilarly for the mean value of the ISE we get 



MISE = 4r < 
M 



p(x)|F(x)-F tr „ e (x)| 2 dx> 



(5) 



Values of e minimising the MISE or the MASE should give 
the optimal representation of the force of the system. 

In the above C is a multiplicative constant, introduced 
to permit comparisons between different mass distributions. 
If we only want to assess the effect of the number of par- 
ticles, as in M96 and sections ^| and |(| we can simply use 
C — 1. On the other hand if more than one configuration 
are to be compared then it may be necessary to use the 
multiplicative constant C. To make this clearer let us con- 
sider two Plummer spheres of the same mass, and of which 
the second one has double the scale length of the first one. 
In that case the optimal softening for the second should be 
double the optimal softening for the first one. By rescaling 
the second Plummer sphere, i.e. multiplying all distances 
by 0.5, we would get for both the same softening, which is 
obvious since they are in fact the same object. It is thus 
preferable to compare softenings after the objects have been 
rescaled to the same mass and size. In a similar way if we 
want to compare the appropriate softening for a Plummer 
sphere and that of a Dehnen sphere we first have to rescale 
one of the two so that the two objects correspond to the 
same mass and size, otherwise the comparison would not be 
meaningful, and would only tell us that bigger objects need 
bigger softenings. Once the objects are rescaled to the same 
mass and size, then the comparison of the optimal soften- 
ings could tell us something about the effect of the central 
concentration. Such a rescaling is, however, not possible in 
all cases. If in a simulation we want to see e.g. the evolution 
of a given Plummer sphere in presence of a given Dehnen 
sphere, then we can not rescale each component separately, 
without changing the problem we are considering. In this 
case we have to search for the softening that would help 
best represent the entire configuration, and then compare 
it with the softening that represents best each of the two 
unsealed objects individually. The above is of course only 
common sense. We have nevertheless explained it in some 
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detail here since otherwise it may be unclear to the reader 
why, in the following sections, in some cases we use a scaling 
and in others we do not. 

We will use the following definition for C which makes 
the MISE and MASE quantities dimensionless: 

C = F' 2 (TZ) 

or, equivalently, 

C = tz a m- 2 

Here M is the total mass in the configuration, 1Z is a char- 
acteristic radius, e.g. the half mass radius, and F(1Z) is the 
force at that radius. Similarly, in such cases, the softening 
should also be scaled with the adopted characteristic radius, 
since it has units of length. 

It is debatable whether in the definition of ISE it would 
have been better to use the numerical or the analytical den- 
sity of the configuration. We have chosen the latter for two 
reasons. One is continuity with the definition of M96, and 
the other is that the calculation of the density of an TV-body 
configuration is not straightforward and may, by itself, in- 
troduce further uncertainties and complications. As will be 
shown in the next section, for sufficiently large number of 
realisations, the present definitions of MISE and MASE 
give the same results and this legitimises our choice for the 
density. 

In all the examples in this paper, unless otherwise 
stated, we have used 6 xl0 6 /TV realisations. Because of the 
high CPU speed of our two GRAPE systems we have used, 
whenever possible, MASE calculations on GRAPE. The 
two notable exceptions are the calculations with alterna- 
tive softening (section uv and the calculations with the stan- 
dard treecode (section g), both of which can not be done on 
GRAPE, and for which we have confined ourselves to radial 
MISE calculations. 

For the treecode calculations on GRAPE we have used 
the software described in detail in A+98. For the stan- 
dard treecode calculated on a workstation we have used the 
Barnes version of the Barnes-Hut algorithm (Barnes and 
Hut 1986) included in the NEMO package (cf. Teuben 1995). 
This version of the treecode can use up to quadrupole terms. 

For the integration along the radius in the radial ISE 
calculation we use the alternative extended Simpson's rule 
(Press et al. 1988), which has an accuracy of 0(N~ 4 ), and 
100 points along the line of integration. For the calculation 
of the radial ISE values for the Plummer sphere we have 
used as an upper limit of the integration L — 20a p , where 
a p is the scale-length of the Plummer sphere. This radius 
contains more than 99% of the mass of the Plummer model, 
while the density has fallen to roughly 3 x 10~ 7 of its central 
value. 



3 OPTIMAL SMOOTHING FOR A PLUMMER 
SPHERE DENSITY DISTRIBUTION 




10 2 10 1 1 10 1 



Figure 1. MASE as a function of the softening e for a Plum- 
mer sphere. From top to bottom the curves correspond to TV = 
30, 100, 300, 1 000, 3 000, 10 000, 30 000, 100 000, 300 000, where TV 
is the number of particles in the realisation of a Plummer sphere. 
The number of realisations taken for the mean is 6 X 10 s /TV. The 
position of the minimum error along a line corresponding to a 
given TV is marked by an X , and the corresponding e value is the 
optimal softening e op t for this number of particles. 

3.1 Dependence of the error on the softening and 
the number of particles 

Following M96 and A+98 we will use, as a first model, a 
truncated Plummer sphere 

p(r) = < 4 ™p v ' p ' 

\ r> R 

where Mt is the mass of the Plummer sphere had it ex- 
tended to infinity, ap is its scale-length and R is its trunca- 
tion radius. For the remainder of this paper, unless otherwise 
noted, we will adopt, without loss of generality, ap = 1, as a 
truncation radius the radius containing 0.999 Mt, and the 
mass within the truncation radius, M(< R) = 1. The corre- 
sponding radial MISE and MASE values were calculated 
in M96, for TV between 30 and 30 000, and in A+98 for TV be- 
tween 30 and 300 000 with direct summation and 1 000 000 
using a tree code. It was found that, as expected, the er- 
ror for a given number of particles TV, be it radial MISE or 
MASE, shows a minimum for a given value of e. In Figure [l] 
we show a number of examples for various values of TV. The 
information is essentially the same as in Figure 4 of A+98, 
except now for MASE rather than radial MISE. 

For small values of the softening the noise dominates 
the error. For this reason the MASE, for such values of the 
softening, decreases steeply with TV, the number of particles 
in the configuration. Conversely for large values of the soft- 
ening it is the bias that dominates. For a sufficiently large 
value of the softening the MASE does not show any depen- 
dence, either on the softening, or on the number of particles. 
In this region the inter-particle forces go as re~ 3 , and tend 
to zero as e — > oo. Thus, for sufficiently large values of the 
softening, the F true (x) term dominates in the difference in 
equations (|), (|), @ and (§) and MISE and MASE tend 
to 
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< / p(x)|F trtle (x)| 2 dx > (6) 
and 

N 

<^Ei Fi ™'Wi 2> ( ? ) 

respectively. These quantities depend only on the mass dis- 
tribution in the configuration, and are independent of both 
the softening and the number of particles, as borne out in 
Figure |l| 

In between the region dominated by the noise and that 
dominated by the bias there is a minimum of the error. Since 
the noise decreases considerably with N, while the bias is 
not affected by it, the minimum error should decrease with 
increasing N and should move to smaller values of the soft- 
ening, as is indeed shown to be the case in Figure |l| 

3.2 The optimal value of the softening and the 
corresponding value of the error 

Let us call e op t the value of e which gives the best repre- 
sentation of the force, i.e. gives the lowest value of MISE 
or MASE (which we will denote hereafter by MISE opt or 
MAS E opt ). Using least square fits, M96 showed that power 
laws give good fits for the values of e op t and MISE op t as a 
function of N, the number of particles in the configuration. 
For N between 30 and 300 000, and C = 1, A+98 obtained 



or of the order of 1.5), where the bias predominates and 
the contribution from the noise is very small. This advan- 
tage is linked to a corresponding disadvantage, namely that 
the CPU time necessary for calculating an ASE is larger 
than the corresponding time for an ISE calculation by a 
factor (JV — 1) /riint, where N the number of particles in the 
representation and n% n t the number of points at which the 
integrand is calculated. This makes the MASE calculations 
prohibitively expensive, even on powerful workstations. On 
the other hand GRAPE-3 and GRAPE-4 are well adapted 
for such calculations. Thus a MASE calculation for 100 000 
particles goes roughly 2000 times faster on our GRAPE-3AF 
system than on its front end, an Ultra 2/200. 

3.4 Comparing results with GRAPE-3 and 
GRAPE-4 

Comparing the MASE results calculated with GRAPE-3 
and GRAPE-4, respectively, we find that they agree very 
well. This could, at first sight, seem at odds with the fact 
that GRAPE-3 has a much more limited accuracy than 
GRAPE-4 that has near 64-bit precision. As argued, how- 
ever, in A+98, the error in GRAPE-3 comes from round-off 
and can therefore be considered as random. Thus when one 
adds the forces from many particles this error cancels out 
and one obtains an accuracy similar to what is obtained 
on GRAPE-4. For this reason we have used both GRAPE 
systems for the calculations given in the next few sections. 



e ovt = 0.987V- U '" b 

and 

MISE op t = 0.22/V"- - 68 

When few particles are used (M96) the exponent of the 
TV dependence of e op t takes a value around -0.28 and that 
of MISEopt around -0.66. If we consider N values between 
10 000 and 300 000 then these values change to -0.23 and 
-0.76, while the asymptotic values for N — > oo are -0.2 and 
-0.8 (A+98). This argues that the above equations are only 
approximate and that in fact the exponents are functions of 
the number of particles considered. Nevertheless, within the 
range of particle numbers used in most collisionless TV-body 
simulations they constitute a very good approximation. 

3.3 Comparing MISE and MASE 

ISE and ASE are, of course, only different ways of calcu- 
lating the same integral. In the case of radial ISE, however, 
the hypothesis is implicitly made that the configuration is 
spherically symmetric, which is obviously true for the con- 
tinuum distribution, or in the limit of an infinite number of 
particles, but not, strictly speaking, in the case of a represen- 
tation of a Plummer sphere with a finite number of particles. 
A different way of seeing the same effect is to say that many 
more points are used in the sampling of the forces in case 
of an ASE than in the case of this ISE. Thus the ASE 
(MASE) errors are smaller than the corresponding radial 
75*75 (MISE) ones for the range of values of e considered 
e.g. in Figure oj, except for the largest values (e greater than 



4 THE CASE OF TWO PLUMMER SPHERES 

In most simulations of galaxies or systems of galaxies it 
is necessary to represent more than one component, each 
having considerably different properties. In order to assess 
the effect of this on the choice of the softening we will in 
this section consider the case of two non-truncated Plum- 
mer spheres. The first one has a scale-length ai = 1 and the 
second one a scale-length a,2 = 0.1. We consider mass ratios 
of the two components Mi /(Mi + M 2 ) = 0, 0.1, 0.25, 0.5, 
0.75, 0.9 and 1. The total mass is in all cases equal to 1, 
the total number of particles equal to 100 000 and the to- 
tal number of realisations equal to 40. We also consider two 
spatial configurations, one in which the two spheres are con- 
centric and the other in which their centers are at a distance 
of 10 length units. For some mass ratios the former configu- 
ration can be considered as representing crudely a halo and 
bulge system, and the latter a halo of a target galaxy with 
a spherical companion. We will discuss further only the case 
of the two concentric spheres, since the two configurations 
give essentially the same results for MASE. As will be later 
shown in the next section it is the distance to the nearest 
neighbours that mainly determines MASE and the corre- 
sponding e op t, and this depends essentially on the density of 
the Plummer sphere with the smallest scale-length, rather 
than on its location. 

Let us first discuss the case of a single simulation in 
which more that one component is present. Since in this 
case we can not apply a separate scaling to each component 
(cf. section |^) we have to take C = 1. Fig ^ compares the 
MASE curves for the seven mass ratios under considera- 
tion, and shows that more concentrated configurations have 
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Figure 2. MASE as a function of e for 100 000 particle config- 
urations. The lower heavy line corresponds to a Plummer sphere 
of unit scale-length and the upper one to a Plummer sphere of 
0.1 scale- length. The remaining ones correspond to cases where 
both Plummer spheres are present, with mass ratios as given in 
the figure. The X symbols show the positions of the minimum 
on each curve. As discussed in the text, for these plots we have 
used C = 1, since we are considering all components present in 
the same simulation. 
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Figure 3. Optimal softening (open circles) and corresponding 
MASE pt (x) as a function of the number of particles in the 
least dense component. The points represented here are obtained 
from calculating the minima of the curves of the previous figure. 



larger errors and smaller corresponding e opt than less con- 
centrated ones. The difference between the respective curves 
is considerable, since a change of a factor of 10 in the scale- 
length makes a change of more than 10 4 in MASE op t and of 
roughly 10 in e opt . This figure also shows that in an iV-body 
simulation with many components the force on the densest 
components will be less well represented than the force on 



Figure 4. Optimal softening (open circles) and corresponding 
value of the MASE ( X ) as a function of the number of particles 
in the least dense component. 



the more extended ones. The total error is bigger when the 
percentage of mass in the denser component is larger. We 
also note that the increase in the error obtained by sub- 
stituting 10% of the particles in a loose configuration with 
a denser one is very big, whereas the decrease in the er- 
ror obtained by substituting 10% of the particles in a dense 
configuration by a looser one is quite small. This is further 
stressed in Figure |^, which is obtained from the minima of 
the curves in Figure ^[ and shows how the optimal soften- 
ing depends on the percentage of mass in the least dense 
component. Similar calculations (not illustrated here) show 
that there is hardly any difference between the cases where 
the two Plummer spheres are concentric and the case where 
they are separated, which leads to the conclusion that the 
densest part influences the result always in the same way, 
independent of its location in the configuration. 

Now let us consider a different question and compare 
seven different simulations, in each of which one of the above 
seven configurations is present. Since now the scaling can be 
applied independently to each of the configurations, in or- 
der to compare these cases between them we need to use 
the weighted versions of the MASE definition and soften- 
ing. Since the total mass of all the configurations is the same, 
the only factor that is changing from one configuration to 
another is the half mass radius, which takes respectively the 
values 0.13, 0.14, 0.18, 0.45, 0.97, 1.18 and 1.30. Now the re- 
sult of the least dense (op = 1) and most dense (ap = 0.1) 
configuration are of course identical. This simply reflects the 
fact that the densest Plummer sphere can be represented as 
well as the least dense one, provided one uses appropriately 
weighted softening values. This is clear also also from Fig- 
ure H which shows MASE opt and function of the 
percentage of mass in the least dense component. This fig- 
ure also shows that, from the configurations analysed here, 
the largest error corresponds to the case with 25% of the 
particles in the more concentrated Plummer sphere. This is 
also the configuration for which e opt is minimum. 
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Figure 5. Comparing the density distributions of the Plummer 
sphere (solid line) and the Dehnen 7 = sphere (dot-dashed line) 
used in section 



Figure 7. Comparing the weighted MASE as a function of e 
for six hundred 10 000 particle representations of a truncated ho- 
mogeneous sphere (dotted line), a Plummer sphere (dashed line) 
and a Dehnen sphere of index 7 = (solid line). 
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Figure 6. Comparing the unweighted MASE as a function of e 
for six hundred 10 000 particle representations of a truncated ho- 
mogeneous sphere (dotted line), a Plummer sphere (dashed line) 
and a Dehnen sphere of index 7 = (solid line). 
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Figure 8. MASE op t as a function of N for a truncated homoge- 
neous sphere (dotted line), a truncated Plummer sphere (dashed 
line) and a truncated Dehnen sphere of index 7 = (solid line). 




5 THE EFFECT OF CENTRAL 
CONCENTRATION 

5.1 Comparison of three different spherical 
distributions 

We have so far considered the case of a Plummer 
sphere, a mass distribution frequently used in astrophysics. 
In this section we will consider two other density distribu- 
tions, a truncated homogeneous sphere and a Dehnen sphere 
(Dehnen 1993). The former is less centrally concentrated 
than the Plummer sphere and the latter more, so that by 
comparing the respective MISE or MASE we can test the 
effect of central concentration on the optimal e and on the 



corresponding accuracy. Since both Plummer and Dehnen 
spheres extend to infinity, we have introduced in both cases 
a cut-off radius, taken so that the mass within that radius 
is equal to 0.999 times the total mass. 

The density profile of the homogeneous sphere is 



p(r) 



3M T /4nR 3 r <R 
r> R 



where R is its outer radius and Mt its total mass. For the 
Dehnen sphere we have 

(3- 7 )M T 

p(r) = i 47r rT(r + a D )(4-T) 







r < R 
r>R 
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Figure 9. e op t as a function of N for a truncated homogeneous 
sphere (dotted line), a truncated Plummer sphere (dashed line) 
and a truncated Dehnen sphere of index 7 = (solid line). 



Figure 11. Optimal softening as a function of mean distance to 
the 6th nearest neighbour. Values for the homogeneous sphere are 
given by open circles and a dashed line, values for the Plummer 
sphere by asterisks and a full line and values the Dehnen sphere 
by X symbols and a dot-dashed line. 
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Figure 10. Histogram for inter-particle distances for the homo- 
geneous sphere (dotted line), the Plummer sphere (dashed line) 
and the Dehnen 7 = sphere (solid line). 

where Mr is its total mass, au is its scale-length and 7 is 
its concentration index, determining the slope of the density 
at the origin. For the examples discussed below we took 
an =0.1 and 7 = 0. 

The truncation radii, containing 0.999 of the total mass, 
are equal to 299.8 and 38.71 for the Dehnen and Plummer 
configuration respectively. For the homogeneous sphere we 
have also taken R — 38.71 and in all three models we have 
taken the mass within the truncation radius to be equal 
to 1. The Plummer and Dehnen density distributions are 
compared in Figure ^[ 

Figure [] and (?) compare the results for MASE obtained 
from six hundred representations of 10 000 particles each. 
In the former we use the un- weighted definition of MASE 



(i.e. with C — 1) and of the softening, and in the latter the 
weighted definition. We note that, in both cases, the Dehnen 
sphere, which is the most concentrated of the three configu- 
rations, requires smaller softening values and is always less 
accurately represented than the two more spread out con- 
figurations, in good agreement with the results of section ^j. 
The differences, however, are much more important in the 
comparison of the un-weighted MASE. Furthermore in this 
case the differences between the Plummer and the homoge- 
neous sphere are very large, while when the comparison is 
between the weighted functions the difference is very small. 
Again this is in agreement with the results of section ^. If 
we simulate all three spheres in the same configuration we 
have to use for all cases the same softening and C = 1. Now 
the forces of the Dehnen sphere will be very badly repre- 
sented and those of the homogeneous sphere very well. The 
situation is totally different if we are simulating one of those 
spheres only, because then we have to calibrate the lengths 
appropriately, so that all systems have the same half-mass 
radius. Now the Dehnen sphere will do much better than 
in the un-weighted case, and the homogeneous sphere much 
worse. Thus the weighted MASE opt is respectively 0.001, 
0.001 and 0.01 for the homogeneous, Plummer and Dehnen 
spheres. In particular the differences between the Plummer 
and the homogeneous sphere are very small. The MASE op t 
for the Dehnen sphere is considerably larger than for the 
other two, presumably because, in the Dehnen sphere, we 
are trying with a single value of the softening to accommo- 
date both very dense and very sparse regions. 

Figures |^ and ^ compare weighted MASE opt and e opt 
as a function of N for the three configurations and confirm 
the trend seen in Figure [?]. The dependence of MASE op t 
and e op t on N for all three density distributions can be rep- 
resented by power laws 

MASE opt = BN b (8) 
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Table 1 . Coefficients of the power laws in equations and 
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-0.69 
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Plummer sphere 


0.64 


-0.25 


0.94 


-0.72 


(weighted) 










Dehnen 7 = sphere 


0.32 


-0.27 


7.5 


-0.69 



(weighted) 

and 

e opt = AN a (9) 

The values of the coefficients are given in Table [l| Since 
the exponent in the above dependences depends somewhat 
on the range of TV used (cf. A+98 and section [3.2] ) , we used 
for the three mass distributions the same values of TV for the 
linear fits, namely the values for TV = 1 000, 3 000, 10 000, 
30 000, 100 000 and 300 000. 

The table and figures show clearly that more centrally 
concentrated configurations need smaller values of the soft- 
ening for an optimal representation of the force, and the 
precision achieved is never as good as for less centrally con- 
centrated configurations. The minimum error, MASE op t, 
decreases somewhat faster with TV for the case of the Plum- 
mer sphere, but the differences are small. 

Figure |l^ compares the histograms of the inter-particle 
distances for the three models after they have been rescaled 
so that the half-mass radii are in all three cases the same. 
This scaling is appropriate for understanding the results ob- 
tained with weights. Each histogram has been obtained from 
ten 10 000 particle realisations of each model. We note that 
the peak of the histogram is nearest to the center for the 
Dehnen sphere, followed by the Plummer sphere, while the 
peak of the histogram for the homogeneous sphere is yet 
further out. It is thus expected that there are more particles 
very close to each other in the Dehnen sphere than in the 
Plummer one, and even more compared to the homogeneous 
one. However the Dehnen sphere has also more particles with 
very large inter-particle distances, while these are fewer for 
the Plummer sphere and even more so for the homogeneous 
sphere. 

5.2 Extension to other configurations 

The results shown in Figure ^| give us the optimal value 
of the softening, which gives the best representation of the 
forces, as a function of the number of particles TV. However 
they also show that a value of softening which is optimal 
for one type of mass distribution is not necessarily optimal 
for another, and that the optimal value depends strongly on 
the central concentration of the distribution. Thus, in order 
to find e op t for a mass distribution other than the three dis- 
cussed here, one can either do the full calculations of the 
MASE, as above, or use the above results to obtain rough 
estimates. Since the former is rather demanding, we would 
like, in the remaining of this section, to see whether it is 



possible to obtain some, albeit crude, estimates of the opti- 
mal softening as a function of quantities linked to the mass 
distribution, but more straightforward to calculate than the 
MASE. 

As we have already seen, smaller values of the softening 
are necessary for more compact configurations or for rep- 
resentations with a larger number of particles, i.e. in cases 
where the particles are closer together. This suggests that 
some measure of inter-particle distances could be used for 
estimating the optimal value of the softening. Furthermore, 
since MASE is more dependent on the accuracy of the forces 
to nearby particles, we will try using the distances of the few 
nearest neighbours. Let us thus, for a given configuration, 
measure, for every particle, the distance to its twelve nearest 
neighbours. Now we need to average this over all particles 
in the configuration, in order to obtain, for the whole con- 
figuration, the mean distance to the nearest neighbour, the 
mean distance to the second nearest neighbour etc., up to 
the mean distance to the 12th nearest neighbour. A stan- 
dard arithmetic average would not be appropriate for this. 
This can be understood if we mentally add to the configura- 
tion a single particle, located so far from it that it can, for 
argument's sake, be considered at infinity. This new particle 
will not influence the value of MASE, the value of e opt that 
should be used, or the accuracy in the force calculations. 
On the other hand it will influence the mean distance to the 
fcth neighbour. It is thus not reasonable to expect a close 
relation between e opt and the arithmetic mean of the fcth 
closest neighbour of all particles. For this reason, instead 
of the standard arithmetic mean, we will prefer using the 
harmonic mean 

JV 

rk,meanl = (TV" 1 r^j)" 1 
i=l 

or, since the force is inversely proportional to the square of 
the distance, 

N 

r k ,mean2 = (TV" 1 ^ r^)' 1/2 
i=l 

where the summations are over all the particles in the con- 
figuration and fc = (1, 12) . In this way we obtain for a 
given configuration the mean distance to the nearest neigh- 
bour, the mean distance to the second nearest neighbour 
etc., up to the mean distance to the 12th nearest neigh- 
bour. It is possible to diminish the noise in such calculations 
by considering several realisations of the same configuration 
and then mean r k , mean i or r k , m ean2, now using simple arith- 
metic means over all realisations. We calculated these mean 
inter-particle distances as a function of TV for the three den- 
sity distributions considered above, using 3 xl0 6 /TV reali- 
sations. In all cases the dependence is roughly linear on a 
log-log plane, and thus can be represented by power laws of 
the type 

where r k>mean is equal to r kimeanl or r k , mean2 . In both cases 
the value of a is around -0.33 or -0.34 and does not depend 
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Table 2. Coefficients of the power laws in equation ( |lo| ) 
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0.95 


0.78 


0.55 


0.76 


0.31 


0.83 


3 


0.59 


0.78 


0.35 


0.76 


0.19 


0.83 


5 


0.50 


0.78 


0.30 


0.76 


0.16 


0.83 


7 


0.45 


0.78 


0.28 


0.76 


0.15 


0.83 


9 


0.41 


0.77 


0.26 


0.76 


0.14 


0.84 


11 


0.39 


0.77 


0.25 


0.76 


0.13 


0.84 



much on which of the twelve nearest neighbours is consid- 
ered. The latter is not true if we use standard means. The 
values of A depend of course on which of the nearest neigh- 
bours is chosen. 

We have thus obtained so far, for a given density dis- 
tribution and a given number of particles N, the average 
distance to the fcth nearest neighbour, for k between 1 and 
12, as well as an optimal softening e opt (cf. Figure^). From 
these two we can eliminate the dependence on and obtain 
the dependence of e opt directly on the average distance to 
the fcth nearest neighbour. This is given in Figure [n] for the 
sixth neighbour, after both distances and the softening have 
been weighted appropriately, as discussed in section Pj. We 
have repeated this exercise for all other values of k between 1 
and 12, but since the results are similar we do not reproduce 
them here. Again a power law gives a good representation 
of the dependences 

The values of A and a, for various neighbours are given in 
Table |. 

Figure [Tl] can give some, albeit rough, estimate of the 
optimal softening to be used, once an rk,mean has been calcu- 
lated, since the differences between the three models are less, 
or of the order of, 0.5 dex. It is, however, possible to narrow 
down the prediction further. Indeed the three dependences 
are ordered as a function of the central concentration of the 
corresponding models, less concentrated models correspond- 
ing to higher softening than more concentrated ones, in good 
agreement with what was previously discussed in this sec- 
tion. Thus comparing the central concentration of the new 
model, whose optimal softening we want to estimate, with 
that of the three considered here, should probably narrow 
the estimate to 0.2 or 0.3 dex. 
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Figure 12. MASE as a function of softening for a homogeneous 
Ferrers ellipsoid of axial ratio 10:10:1. The various curves corre- 
spond to a different number N of particles in the configuration; 
10 000 (solid line), 30 000 (dashed line), 100 000 (dot-dashed line), 
and 300 000 (dotted line), respectively. 




e 



6 THE EFFECT OF TRIAXIALITY 

So far we have considered only spherical objects. By 
using, however, the MASE rather than the radial MISE 
accuracy estimator it is possible to consider non-spherical 
configurations. As an example we will in this section consider 
a Ferrers' ellipsoid (Ferrers 1877), a distribution often used 
for modeling bars. The corresponding density is 




Po(l-g 2 ) n if 9 <l 



where g 2 = x 2 /a 2 + y 2 /b 2 + z 2 /c 2 , a > b > c are the major, 
intermediate and minor axes, po is the central density of the 
ellipsoid and n is an index determining the mass distribution 
in the ellipsoid. We calculated MASE values for four such 



Figure 13. MASE as a function of softening for four Ferrers' 
ellipsoids of axial ratio 1:1:1 (solid line), 12:4:3 (dashed line), 
10:10:1 (dash-dotted line) and 74:74:1 (dotted line). In all cases 
the number of particles is equal to 100 000 and the number of 
realisations is equal to 60. 

ellipsoids with n — and a : b : c — 1:1:1, 12:4:3, 10:10:1 and 
74:74:1 respectively. These represent objects common in an 
astrophysical context: a sphere, a bar, a thin and a very thin 
disc. Undoubtedly the disc component in real galaxies is not 
as thin as the thinnest of the two discs that we are using here, 
but we have added on purpose this rather extreme ellipsoid 
to see how an extreme thinness would affect the values of 
MASE. To facilitate comparisons, the four ellipsoids have 
been taken to have the same volume and mass (equal to 
666.666 and 1 respectively), so we can take C = 1. 
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Figure 14. Histogram of inter-particle distances of all pairs of 
particles in Ferrers' homogeneous ellipsoids with axial ratio 1:1:1 
(dot-dashed line), 12:4:3 (dotted line), 10:10:1 (dashed line) and 
74:74:1 (solid line). Each histogram was obtained from ten 10 000 
particle realisations. 
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Figure 15. Bias for a softening tending to infinity as a function of 
the axial ratio c/a of the ellipsoid. Asterisks correspond to oblate 
ellipsoids and open squares to prolate ones. Open circles and X 
correspond to axial ratios a : a/2 : c and a : a/4 : c respectively. 
All ellipsoids have the same volume as the four studied in more 
detail in this section. 



Figure |l2| shows MASE as a function of the softening 
for a Ferrers homogeneous ellipsoid of axial ratio 10:10:1 and 
different values of N. The general outline of t he c urves is the 
same as for the spherical objects (cf. section |3.l| ). As is the 
case for the spherical distributions e opt can be represented 
by a power law dependence on the number of particles iV. 

In order to evaluate better the effect of triaxiality on 
the error in the force calculation we plot in Figure Il3l the 
MASE as a function of e for the four ellipsoids under con- 



sideration. In order not to burden the figure we plot only 
one value of N, in this case 100 000. Similar results have 
been obtained for the other values of N considered. Let us 
first concider the flat part of the curve, where the MASE is 
practically independent of the softening, and which occurs 
for large values of e. For brievity we will hereafter call this 
the bias part of the curve. Two effects in particular are clear 
from Figure |l^ about the bias part. The first concerns for 
how big a value of the softening this part is attained and 
the second what the value of MASE on this section is. 

Figure [l^ shows that the bias part occurs for relatively 
smaller e values when the shape is spherical, and consid- 
erably larger ones when the departure from sphericity is 
large. This can be understood with the help of Figure 
which compares histograms of the inter-particle distances 
for the four ellipsoids under consideration. We note that, 
as expected, the higher the departure from sphericity, the 
larger the percenta ge o f large inter-particle distances. As 
discussed in section 3.1, the bias dominates when the soft- 
ening is larger than most inter-particle distances. According 
to Figure |l4l this is expected to happen for larger soften- 
ings for objects that depart more from sphericity and this is 
indeed verified in Figure |l3[ 

The second effect that is clear from Figure |l^ is that the 
value of the bias, or, in other words, the values of MASE 
on the bias part of the curve, also depends on the depar- 
ture from sphericity. In particular the larger the departure 
from sphericity, the smaller this value is. Th is c an again 

We can 



5.1 



be understood with the discussion in section 
simply calculate the value of the bias at infinite softening 
from equations Q or (Q) by integrating the force of a ho- 
mogeneous prolate spheroid appropriately over its volume. 
We find the highest biases for the most spherical objects, 
and the values decrease as the departure from sphericity in- 
creases, in good agreement with the results of Figure |l3| 
Some examples are shown in Figure |l5|. We see that both 
oblate and prolate objects have high values of the bias when 
they are nearly spherical. When they are far from spherical 
then the bias for the prolate cases is somewhat higher than 
the bias for the corresponding oblate ones. Figures [l3] and 
[lfj show that the effect of shape on the bias can be quite 
important, changing the corresponding values of MASE by 
over an order of magnitude. 

Now let us turn to the values of MASE corresponding 
to smaller softenings. Figure |l3| shows that they do not de- 
pend much on the shape. This can be understood since, for 
such values of the softening, it is the noise that affects the 
error most, and this should not depend on the shape of the 
object. In particular for values around e op t this dependance 
is negligible. Thus we can conclude that the shape of the ho- 
mogeneous Ferrers ellipsoid does not influence much either 
the value of e opt which brings the best representation of the 
forces, or the corresponding value of the error MASE opt . 



7 DIFFERENT FORMS OF SOFTENING 

In all calculations presented so far we have used the 
standard Plummer softening, introduced in section ^. There 
are, however, many alternative ways of introducing the 
smoothing. For example, instead of using the second power 
of the softening and of the inter-particle distance in equation 



© 1999 RAS, MNRAS 000, 0-|l| 



12 



0.8 - 



0.6 



0.4 - 



0.2 




10 



10 



10 



;:. r.-- e ..o .<>..,,, .a. .,, ^ . 



4 6 
softening power p 



Figure 16. Comparison of the forces calculated with three dif- 
ferent types of softening to the Newtonian force (solid line). The 
Plummer softening is given by a dashed line, the power softening 
with p=4 with a dot-dashed one, and the spline softening with a 
dotted line. 



Figure 18. As in the previous figure, but for the minimum value 
of the radial MISE. 
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Figure 17. Optimal softening as a function of the exponent p 
in the power law softening. The solid line and stars correspond 
to N = 10 000, the dashed line and open circles correspond N = 
30 000, and the dot-dashed line and X symbols to N = 100 000 
(cf. text). The corresponding values for the spline softening arc 
given by horizontal lines of the same type. 



(^), one could use any other value, including non-integer. 
Another alternative is to use, instead of the force in a Plum- 
mer sphere, the force in a sphere of constant density (e.g. 
Pfenniger & Friedli 1993, Palous, Jungwiert & Kopecky 
1993), or of any other radial dependence (e.g. Palous, Jung- 
wiert & Kopecky 1993). Alternatively one could substitute 
the force between two point masses for the force between 
two spheres. A few density profiles of the mass within the 
sphere lend themselves to an easy calculation of the potential 
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Figure 19. Optimal softening and corresponding minimum value 
of the radial MISE as a function of N for a Plummer sphere. 
This figure compares the results for a spline softening (solid line), 
a power softening for p = 2 (dashed line), and a power softening 
for p = 5 (dot-dashed line). 



and corresponding force, in particular certain radial profiles 
given by polynomials. Such examples have been discussed 
by Hockney & Eastwood (1981) and by Dyer & Ip (1993). 

We will discuss in this section some alternative types 
of softening and compare how well they fare in the MISE 
test. Since both GRAPE-3 and GRAPE-4 calculate the force 
only with a Plummer softening, we performed all calcula- 
tions in this section on workstations, using radial MISE, 
rather than MASE, and 64 bit precision. The integration 
in radial MISE was done as described in section ^. 

As a first example of an alternative force calculation 
we will use the spline approximations given by Hernquist & 
Katz (1989), namely F = — mr/(r), where 
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l/e J [(4/3) - (6/5)^ + (1/2) U J ] 
f{r) = ■> l/r 3 [-l/15 + (S/3)u 3 - 3u 4 + (6/5)u 5 
1 I/- 3 



< u < 1 

■ (i/e)u 6 ] i < u < 2 

u > 2 



Table 3 . Coefficients of the power law fits 



and it = r/e. 

Note that for r > 2e the value of the force in this ap- 
proximation is exactly the Newtonian force. In other words 
contrary to the Plummer softening, and to the power-law 
softening introduced below, the spline softening is approxi- 
mate only for small distances. 

As a second example we will consider an extension of 
the Plummer softening to values of the exponent other than 
two. This can be given by 



Fi = G 



N i 

J'=l 



- Xi)|Xi 



(|x 1 -x j |p + £P) 1 /p+i 



(11) 



which for p=2 gives back equation (|l|). A point to note is that 
for all values of p, including the commonly used value p = 2, 
these forces tend to the Newtonian one only asymptotically, 
i.e. even at large distances there is a finite, albeit small, 
difference between the results they give and those of the 
Newtonian force. In other words they introduce a small but 
non-zero smoothing even at large distances, where it isn't 
necessary. 

Figure ^6| compares the amplitude of the non-softened 
Newtonian force with those of the Plummer, the p = 4 
power-law and spline softened forces. For this figure we have 
taken the softening as well as the masses equal to unity. We 
note that the force that approximates best the Newtonian 
one is the spline, followed by the higher power softening, 
while the Plummer softening does the least well of the three. 
Thus the forces agree better than 5% with the Newtonian 
one for distances larger than roughly 1.3, 2.2 and 5 softening 
lengths, correspondingly for the spline, p=4 and Plummer 
softening. 

In order to assess how well each softening can represent 
the forces in a given mass distribution we calculated radial 
MISE values for different values of the softening, for dif- 
ferent number of particles in the configuration, and, in the 
case of the power softening, different values of the exponent 
p, roughly in the range from 2 to 8. For all types of softening 
the MISE as a function of e curves are very similar to those 
of Figure [l], so we will not repeat them here. 

Figure |l7] compares the optimal softening {(.opt) for the 
power law softening - as a function of the exponent p - and 
for the spline softening. The values were obtained from six 
hundred realisations of a Plummer sphere of 10 000 particles 
each, two hundred realisations of 30 000 particles and sixty 
100 000 particle realisations. We note that e ov t increases 
with p. Thus comparing p — 2.0 to p = 7.0 we find an in- 
crease in e opt of roughly a factor of two. The e op t value for 
the spline is somewhat larger than that of the highest p 
values. As expected, the optimal softening decreases with 
increasing number of particles N. The A(loge op t) does not 
depend notably on the power p. 

Figure nM compares the minimum value for radial 
MISE (M I S E opt ) for the same cases as the previous figure. 
We note that MISE op t decreases with p. Thus comparing 
p — 2.0 to p — 7.0 we find a decrease in MISE op t of roughly 
30%. The MISEopt value for the spline is of the order of 
that of the highest p values. As expected, the corresponding 
minimum errors decrease with increasing number of parti- 
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cles N. The A(logMISE op t) does not depend notably on 
the power p. 

Figure |l^ compares the optimal softening and the cor- 
responding radial MISE values as a function of iV for the 
spline softening and for the power softening for exponents p 
— 2 and p = 5. Power laws are satisfactory approximations 
in all cases, given by 



eopt = AN a 



and 



MISEopt = BN 

The values of the coefficients are given in Table |S| The 
small differences between the coefficients for the Plummer 
sphere and Plummer softening given here and those given 
in Table |l| are due to the fact that here we have used radial 
MISE while in Table | MASE. From this table, as well as 
from Figure jH| we see that the MISE op t as a function of 
N is more or less the same for the p = 5 and the spline. The 
p — 2 case gives somewhat bigger values of MISE op t- 

All the above argue that the spline softening as well 
as the higher values of the power in the power softening 
give a better representation of the force than the standard 
Plummer softening. The difference, however, is not as big as 
one could have inferred from Figure [l(| since some of the 
difference is compensated by an adjustment in e ov t- Thus 
Figures |l7| and ^ argue for an improvement of 30%, with 
corresponding changes of e op t of a factor of two. This im- 
provement is nevertheless non-negligible, since it would take 
an increase of the number of particles of roughly 70% to 
achieve it (cf. section [j.2| ). The fact that the corresponding 
value of e op t is higher is also an advantage, since, for equally 
good representations of the forces in the mass distribution, 
a larger softening allows for large time-steps, and therefore 
shorter CPU execution times. Last but not least the spline 
softening necessitates considerably less CPU time per call. 
This gain in time depends on whether one programs in for- 
tran or C, on what the exponent of the power is, and on the 
compiler used. We have found ratios roughly between 2 and 
10. 



8 TREECODE 

By making some modifications to the standard treecode 
(Barnes and Hut 1986) it is possible to implement it on a 
GR APE system (Makino 1991). In particular the tree should 
not be descended for each particle separately, but for blocks 
of particles, as initially proposed by Barnes (1990). Increas- 
ing the number of particles in the block makes the interac- 
tion list longer and the treecode more accurate. The par- 
ticular implementation on the Marseille GRAPE systems 
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is described in A+98, together with some discussion on its 
performance and accuracy. 

We have made calculations of MASE using the GRAPE 
treecode and radial MISE with the standard one, for var- 
ious values of the tolerance and the number of particles 
TV. The differences between the values corresponding to the 
same number of particles and different tolerances is rather 
small. In particular the values obtained with a tolerance of 
0.5 or 0.7 are very near those obtained with the direct sum- 
mation. Only for tolerances larger than 1 do the MASE 
values increase significantly, and even so the differences with 
the direct summation are always considerably smaller than 
those obtained by changing the number of particles by fac- 
tors as those considered e.g. in Figure El. 



9 SUMMARY AND DISCUSSION 

In this paper we have discussed the value of the softening 
that allows us to best approximate the true forces within a 
given mass distribution in an N-body simulation. 

We have first worked with the Plummer sphere and con- 
firmed previous results that, for a given number of particles 
TV, there is an optimal softening which gives the best approx- 
imations to the forces. For smaller values of the softening the 
noise introduces errors, while for larger there is a systematic 
bias from the Newtonian force results. We calculated the de- 
pendence of the optimal softening on the number of particles 
TV and confirmed and extended the results of M96 and A+98. 
We compared the results obtained for integrated and aver- 
age square errors and found that the latter were consider- 
ably less noisy, as could be expected since they involve many 
more samplings. For this reason MASE calculations require 
considerably more CPU time. Since we have at our disposal 
two powerful GRAPE systems we performed MASE calcu- 
lations wherever this was possible, i.e. in all cases with direct 
summation and the standard Plummer softening. We used 
radial MISE calculations for the other softenings. Some of 
the treecode calculations were done on a GRAPE and some 
on a workstation. We also find that results obtained with 
GRAPE-3 agree very well with those obtained on GRAPE- 
4, despite their differences in accuracy. This is due to the 
fact that the errors in GRAPE-3 are due to round-off and 
can thus be considered as random (cf. A+98). 

We then worked with other density distributions and 
found that the density, and the central concentration, have 
a large influence on the optimal softening and the corre- 
sponding errror MASE opt . We first examined the case of 
two concentric Plummer spheres of different scale- length, 
and found that the existence of a dense sphere has a large 
influence on the e opt and MASE opt . We then compared two 
other density distributions to the Plummer one, namely the 
homogeneous and the Dehnen sphere. The former is much 
less centrally concentrated and the latter much more than 
the Plummer sphere. This confirmed the importance of cen- 
tral concentration on the optimal softening. All our results 
show that denser configurations necessitate smaller soften- 
ings and are never as well represented as less dense ones, i.e. 
they have considerably larger values of MASE opt . 

Since the choice of the optimal softening depends on 
the configuration under consideration, and it is rather cum- 
bersome to perform a MASE study for every different con- 



figuration, we propose a simple way of obtaining a rough 
estimate of this optimal value. We show that it depends 
on the mean distance of the nearest neighbours, a quantity 
which is much easier to calculate or estimate than MASE. 
The precision of this rough estimate should be sufficient in 
many cases. 

We next examined the influence of the shape of the 
object, with the help of Ferrers ellipsoids of different ax- 
ial ratios. Although the influence is large on the bias, it 
is very small for the optimal softening and corresponding 
MASE op t. 

We then examined two alternative types of softening. 
One is a power-law softening in which the value of the ex- 
ponent can have values different than 2, and the other is a 
spline softening. Large values of the exponent as well as the 
spline necessitate a larger value of the softening and give a 
smaller value of the MASE opt . Since the higher values of 
the exponent necessitate more CPU time than the spline, 
it is the latter that provides the best ratio of accuracy to 
CPU time. The difference, nevertheless, with the standard 
Plummer softening is not very large. 

The treecode results are somewhat less accurate than 
those obtained with direct summation and the accuracy de- 
creases with increasing tolerance. The differences, however, 
are small. Thus in order to obtain more accurate results 
within a given CPU time it should be preferable to increase 
the number of particles TV, rather than decrease the open- 
ing angle. This is particularly true for the standard treecode, 
where the dependence of the CPU time on the opening angle 
is stronger (Hernquist 1987, A+98). 

The question we have addressed in this paper is ob- 
viously of interest for TV-body simulations. If the adopted 
value of the softening is too small, then the result will be 
very noisy. In TV-body simulations one considers only one 
realisation, thus the effect of noise can be very acute. On 
the other hand if the adopted value of the softening is too 
large, then the simulation will pertain to another object 
than was initially desired. In certain cases this may be of 
little importance, provided the softening is not excessive. 
This is not, however, always the case. For example if we 
want to check the stability of a model obtained with the 
Schwarzschild method (Schwarzschild 1979, 1982), then it is 
important that the forces approximate as closely as possible 
those of the model. 

The size of structures that need to be analysed has 
sometimes been used in simulations to set the value of the 
softening. Our results, however, show that this may not be 
always possible, and should sometimes be accompanied by 
a corresponding increase in the number of particles. Indeed 
if the size of the structures to be analysed are smaller than 
the e op t, then a better resolution can be reached only by 
increasing the number of particles accordingly. 

It is customary in N-body simulations to use a softening 
which is constant both in time and position. Our results, 
however, show that this practice should be questioned. 

It is easy to implement a softening which is a function 
of time. For example in the case we are simulating a col- 
lapse the softening at the initial stages can be considerably 
larger than during the stages of maximum collapse. This 
could also be said for head-on encounters of galaxies, as e.g. 
in the formation of ring galaxies, where the central concen- 
tration rises considerably over a short period of time. Such a 



© 1999 RAS, MNRAS 000, §4l| 



variation of the softening with time would be easy to imple- 
ment both on standard software and on GRAPE and could 
lead to considerably more accuracy. 

More challenging is to consider a softening which a func- 
tion of position, as well as of time, particularly on GRAPE 
systems. For both GRAPE-3 and GRAPE-4 considerable 
modifications of the software are necessary. The results of 
this paper, however, argue that such an effort could be well 
worth-while. Such a work is in progress (Athanassoula & 
Lambert, in preparation). It should also be noted that, when 
changing the softening, corresponding changes in the time 
step should also be considered. 

The fact that the optimal softening is a function of the 
number of particles used in a simulation affects the CPU 
time necessary for a given simulation. Suppose that we dou- 
ble the number of particles N for a given simulation. Then 
the simulation time will be multiplied by 4 if we are using 
direct summation and roughly by 2 if we are using a treecode 
on GRAPE (A+98). However the optimal softening will have 
also decreased by a factor which for a Plummer type concen- 
tration will be of the order of 0.85. This will mean that the 
time step should also be decreased by a similar factor, and 
the total CPU time of the simulation increased accordingly. 
Thus the total CPU time will increase roughly as TV 2,24 for 
direct summation and as TV 1,24 for the treecode. 

The dependence of the optimal softening on the num- 
ber of particles will also affect the relaxation time. Huang, 
Dubinski & Carlberg (1993) derived the following simple for- 
mula for the relaxation time as a function of the softening 
and the number of particles, assuming that the distribution 
of particles is homogeneous. 

N 

T re l ax — 77; / r . , rT croS s (12) 

8ln(R/e) 

In the above T cross is the crossing time and R is some char- 
acteristic radius of the system. Replacing the softening by 
its optimal value, which is a function of the number of par- 
ticles, we find that the relaxation time will grow somewhat 
less rapidly then expected if the softening was not a function 
of N. 

These are the prices to pay for a simulation with bet- 
ter resolution, with reduced particle noise and with more 
accurately calculated forces. 
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